Real-time simulations of nonequilibrium transport in the single-impurity Anderson model 
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One of the main open problems in the field of transport in strongly interacting nanostructures is the under- 
standing of currents beyond the linear response regime. In this work, we consider the single-impurity Anderson 
model and use the adaptive time-dependent density matrix renormalization group (tDMRG) method to compute 
real-time currents out of equilibrium. We first focus on the particle-hole symmetric point where Kondo correla- 
tions are the strongest and then extend the study of the nonequilibrium transport to the mixed-valence regime. 
As a main result, we present accurate data for the current-voltage characteristics of this model. 



I. INTRODUCTION 

Experiments in strongly interacting nanostructures are of- 
ten well described by relatively simple model Hamiltoni- 
ans. In several cases, these models are integrable and can 
be solved exactly, 1 or they can be studied using powerful nu- 
merical methods, such as the numerical renormalization group 
(NRG). 2,3 These techniques have had enormous success in de- 
scribing the equilibrium properties of these models. How- 
ever, the understanding of open systems out of equilibrium 
remains an extremely challenging area of research. While 
recent studies have succeeded in calculating the current- 
voltage (I-V) characteristics of the interacting resonant level 
model j 4 ' 5 ' 6 ' 7,8 the focus of current research has been devoted 
toward the more difficult single-impurity Anderson model in- 
corporating Kondo correlations. For example, experimental 
results for the finite-bias transport properties of a quantum 
dot still await a complete theoretical explanation . 9 ' 1 °' 1 1 ' 1 2 Be- 
sides the fundamental interest in solving this type of problems, 
Kondo physics is an emergent feature in many experiments on 
nanoscopic systems ] 9 ' 10,11 

Our goal here is to calculate the I-V characteristics of the 
single-impurity Anderson model, described by the Hamilto- 
nian H = H dot + iJi cads >I2 

#dot = V s y]n la + U/2'Y] nxqUxa 

Pleads = ~t' ^ (4,<7 C i,l,<T + hx -) 
l=R,L-a 

Nl.r 
n=\\l=R,L:a 

Hdot contains the gate potential V g and the Coulomb repul- 
sion U acting on the dot. The first term in i?i ca ds is the hy- 
bridization of the dot with the conduction band, and the last 
term describes the leads. is the number of sites on the 

left (right) lead, with N = Nr, + Nr + 1. c\, a is a fermionic 
creation operator acting on the dot (v = 1) or a site i in the 
left (right) lead (v = L(R), i), with a spin index a =] , j. 



The Kondo resonance is a signature of the hybridization of 
the leads with a single electron on the dot Electrons from 
the leads screen the extra spin, forming an overall singlet, 
and the corresponding Kondo cloud^ This collective state 
can be suppressed, and eventually destroyed, by temperature, 
magnetic field, or a large bias between the leads. The ulti- 
mate goal is to understand the subtle mechanisms behind the 
suppression of the Kondo resonance, and how this response 
translates into the I-V characteristics of the system. In par- 
ticular, when a finite voltage bias is applied, results for the 
density of states (DOS) at the particle-hole symmetric point 
V g = -U/2 seem to indicat e) 14 - 15 - 16 ' 17 - 18 that the Kondo peak 
smoothly decreases, while the spectral weight is transferred 
to satellite peaks centered at the chemical potentials of the 
leads at frequencies to w ±V/2. At large biases, these peaks 
merge with broad shoulders at u> ~ ±U /2i 14 ' 15 - 16 - 17 ' 18 The 
possible existence of a negative differential conductance in 
the high-voltage regime of nanostructures is a topic of intense 
discussion! 4 - 

Beyond perturbative limits in U or bias voltage V, the 
properties of the steady state, including average currents, 
are still under scrutiny. Therefore, to gauge the validity of 
approximate analytical schemes used to study nonequilib- 
rium transport, it is highly desirable to have accurate nu- 
merical results available as a benchmark. Several efforts to 
work out the nonlinear transport regime of this model have 
been made using perturbativ o 14 ' 19 - 20 as well as renormalization 
schemeS j 21 - 22 ' 23 ' 24 - 25 ' 26 the non-crossing approximationriH 
flow-equations* 2 - 7 - real-time path integral approaches, 28 Quan- 
tum Monte Carlo (QMC) / 5 - 29 ' 30 and the Bethe ansatz^i 
Recently, the time-dependent NRG method has been used 
to simulate aspects of nonequilibrium physics in quantum 
dotsj 16 ' 32 - 33 - 34 In this work, we introduce the tDMRG^£ to 
address two questions: first, what is the I-V characteristics 
at V g — — U/2 and, second, how do nonequilibrium currents 
behave in the mixed-valence regime. Furthermore, tDMRG 
gives access to the transient regime, i.e., information on how 
the steady state is reached, which some of the aforemen- 
tioned methods are, as of now, not designed for. We will 
also provide a technical discussion of our tDMRG calcula- 
tions, in terms of system sizes and accuracy needed to obtain 
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high quality results. We work at values of U/T < 8, where 
r = 2n pi ca & s (E p)t 12 = 2t 12 is the hybridization parameter, 
with p\ aa ds{Ep) the density of states in the tight-binding leads 
at the Fermi energy Ep. Our approach employs a real-space 
representation of the leads with t n = 1, which we find is a 
useful modeling of the leads for bias values V > Tk, where 
Xk is the Kondo temperature. In the case of V g — —U/2, our 
results are in excellent agreement with those from the func- 
tional renormalization group (fRG)22 for U/T < 6. 

The plan of this exposition is the following. In Sec. [II] 
we shall introduce our numerical approach, the tDMRG, and 
we discuss the error sources and the quality of our numeri- 
cal data. In Sec. [Ill] we present our results for the particle 
hole-symmetric point, while in Sec.[lV] we turn to the mixed- 
valence regime. We conclude with a summary and discussion 
of our results in Sec. [V] 



II. NUMERICAL METHOD: DMRG 

The tDMRG method has previously been used in Refs. [37| 
and [3^ to study the linear conductance of the single-impurity 
Anderson model Eq. ((T). In the latter effort, it has been argued 
that a logarithmic discretization^, i.e., t n oc \~ n l 2 is useful 
in capturing Kondo correlations in the linear conductance at 
V g = —U/2 and using this trick, perfect conductance was 
obtained for U/T < 7. Other applications of the tDMRG to 
transport in nanostructures include the (interacting) resonant 
level modelji^ short chains of spinless fermions 4 ^ 4 -!- as well 
as nonequilibrium transport through a point contact in a non- 
Abelian quantum Hall state. 4 ? Moreover, first tDMRG results 
for Eq. ([Til out of equilibrium and at V g = —U/2 were pre- 
sented in Ref.]HbyKirinoef al. for U/T < 2.8 and in Ref. [H 
for U/T = 3.125. We here move substantially beyond the 
scope of these previous studie s 38 ' 43 as we consider larger val- 
ues of U /r, perform a finite-size scaling analysis, and use up 
to m = 2000 DMRG states during the time-evolution (com- 
pared to m = 800 in Ref.l43t). 

In our numerical calculations^!^ we first compute the 
ground state of Eq. ([T]i and then, at time r = + , we apply 
an extra term i^bias and evolve under H + H^i as : 

-Hbias — — n R,n — "J n L,n ■ (2) 

n—1 n—1 

In principle, different spatial forms of V can be usedr^S 
which sometimes affect the transient behavior but leave the 
steady-state currents unchanged. We measure the expectation 
value J(t) of the local current operator acting on the links 
connecting the dot to the leads, averaging over the right and 
left links i^lS We set e = H = 1, All simulations are per- 
formed at an overall half filling of dot and leads. 

In our tDMRG runs, we use a Trotter-Suzuki time evolu- 
tion scheme with time-steps of 0.05 < Srt < 0.15. We find 
a fast convergence with system size on chains with Nl(Nr) 
odd(even)i2i2i Our goal is to achieve a numerical accuracy of 
& J steady /V ~ O.OlGrj, where Go is the conductance quan- 
tum. The steady-state currents are computed by taking suit- 
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FIG. 1: (Color online) Particle-hole symmetric point, V B = —U/2: 
Current J(r)/V vs. time r for U/T = 2 (U = 0M, t' = 0.3535t, 
V s = -17/2) for V/U = 0.4, 0.6, . . . , 2 [JV = 96 (thick lines); 
N = 64 (thin lines); N = 48 (dot-dashed lines)]. The arrows in- 
dicate the time-interval used for the calculation of the steady-state 
current J stC ady = {J( t ))t for V/U = 1.2 from a time average 
( • ) T .— Inset: estimates of the transient time T trans as a function of 
bias for U — t/2 (squares, T = t/4) and U/t — 1 (triangles up, 
T = f/2), both at U/T = 2. 



able time-averages, i.e., J s tcady = (J)t- in principle, the true 
steady state is only reached for N — > oo, while on finite sys- 
tems, the term quasi-steady states may be more appropriate. 
For simplicity, we do not make this distinction here. 

The error sources in tDMRG arei^ 4 ^ (i) truncation errors, 
(ii) the Trotter error, and (iii) the determination of the time in- 
tervals to extract average currents from the real-time data at 
a given bias. Typically, we perform runs with different dis- 
carded weights 8p at fixed U, T, V and N to identify the time 
windows for the averaging. We find (i) and (iii) to be the dom- 
inant sources in the present problem as the Trotter errors can 
effectively be reduced by using small time-steps and this error 
grows only linearly with the simulation time4& The truncation 
errors can be suppressed by decreasing the discarded weight 
Sp (Ref. 1461) during the time evolution. This error source is 
thus well-controlled, too. As for the averaging procedure, 
the determination of steady-state regimes are the least unam- 
biguous for large gate voltages, intermediate biases voltages 

V > 2k, and/or large U/T > 6. There, the accuracy of our 
data for ( J) T /V is worse than 0.01 Go- 

An important aspect of any simulation of time-evolution is 
the entanglement growth (see, e.g., Ref. l47h . This is mea- 
sured through the entanglement entropy Si = — tr(p/ln(p/)), 
where pi is the reduced density matrix for a block of length 
I, as is naturally accessed in DMRG4& As our set-up real- 
izes a so-called global quantum quench, Si can grow signif- 
icantly as a function of time. As a consequence, to keep the 
truncated weight below a given threshold, an increasing num- 
ber of states needs to be kept during the time evolution. For 

V > O.li, it is imperative to work at a fixed Sp (as compared 
to a fixed m) during the time evolution to have firstly, an ef- 
ficient simulation and secondly, a meaningful control over the 
numerical accuracy by studying the scaling in Sp. Typically, 
to meet our target accuracy at V < 2t, a discarded weight of 
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FIG. 2: (Color online) Particle-hole symmetric point, V g = —(7/2: 
(a) Current-voltage curves for U/T = 2, 4, 6 with U — t/2 (squares, 
stars, circles) as well as data for U/T — 2,8 with U = t (triangles 
down and up, respectively). Open symbols: TV = 64, full ones: TV = 
96 sites, shaded ones: TV = 128, stars and triangles down: finite-size 
extrapolations in 1/TV (U/T = 4, 8 only). Dash-dotted line: (J) T = 
2GqV. (b) Enlarged view of the low bias region for U/T = 4 and 
8, {J)r/V vs. V. Solid lines: fRG results from Ref . |H; dotted 
lines: 4th-order perturbation theory results from Ref.0. (from top to 
bottom). For U/T = 8, the arrow in the inset indicates the difference 
between the fRG and tDMRG results on the one hand and those from 
Ref. uA on the other hand. For NRG results for Tk, see, e.g., Fig. 6 
inReTgl 

Sp < 1CP 7 is found to be appropriate. 



III. RESULTS FOR PARTICLE-HOLE SYMMETRY 

After these remarks on the method, we turn to the discus- 
sion of our tDMRG results. The time-dependent currents at 
Vg = -(7/2 are displayed in Fig. Q] for U/T = 2 (U = 0.5t) 
and V/U = 0.4, 0.6, . . . , 2.0. Since we plot J(r) /V, the cur- 
rents decrease as the bias increases. An example for the time- 
interval used to compute the steady-state currents is indicated 
by arrows for the case of V = 1.2U in the figure. While a 
full analysis of the short-time transients and its dependence on 
both V and T (see also Ref. HH) will be presented elsewhere, 
note that the transient times, i.e., the time for the current to 
reach its steady-state plateau, decreases as the bias increases, 
to, e.g., about r ~ 6/t at V > 2U and U = t/2. This is il- 
lustrated in the inset of Fig.Q] where we estimate the transient 
time from the time that it takes the initial maximum in J(r) 
to decay. For this reason, and due to less severe finite-size 
effects, short chains of TV ~ 24, 32, 48 are often sufficient at 
large biases. 

The average currents (J) T are plotted as a function of bias 
in Fig. Ha) for U/T = 2,4,6 (U = 0.5t for all three sets, 
plus additional U = t curves for U/T = 2 and 8). For U/T = 
2 and U = t, the average current increases monotonically 
with bias and eventually decreases. The latter gives rise to 
a negative differential conductance, yet this is simply due to 
the decreasing overlap of the leads' DOS, which have a finite 




N=96 (thin lines) 

1 i i i l i i . i — 

10 20 30 40 

time x t 



FIG. 3: (Color online) Mixed-valence regime: Current J(t)/V vs. 
time t for U/T = 2 (U = t, t' = OM, V g = 0) and bias values 
V/U = 0.3, 0.8, 1, 1.2 [TV = 96 (thin lines), TV = 64 (thick lines)]. 
Arrows indicate the interval used for the averaging for AV = 0.8U 
and TV = 64. The thin solid line for V = 1.2(7 is a fit of a + 
bcos(ujr + cj>) to the corresponding TV = 64 tDMRG data (dotted 
line). 



bandwidth of At in our case. For that reason, namely the finite 
band-width, the current vanishes in the limit of V — > oo in 
our simulations. 

We can estimate at which voltage the band-curvature and 
the finite bandwidth, and thus the deviation from the wide- 
band limit, start to play a role by varying U while keeping 
U/T fixed. Comparing U = 0.5i (squares) and U = t (tri- 
angles up) results in Fig. |2ja) at U/T = 2, we find that the 
smaller Us give the best agreement with fRG 22 (thick, solid 
lines) at large voltages. Note that the fRG calculations at 
these small values of U/T ~ 2 are expected to be very re- 
liable. We thus mostly use U = 0.5i for our runs at larger 
values of U /T as we can then compare to the fRG at larger 
values of V. This has the disadvantage that the transient time 
increases with increasing 1/T (corresponding to a decreasing 
U at a fixed ?7/T;see the inset of Fig.Q~|i, and therefore, for the 
runs at [7/T = 8 where we are interested in the behavior at 
intermediate bias values, we use U = t. 

At U/T = 6, a tendency of separating the narrow Kondo 
peak at small bias values from the Coulomb blockade peak at 
V/U > 0.5 is visible. Unfortunately, at biases V ~ 0.2(7, we 
observe a crossover in the monotonic behavior of our finite- 
size data, with (J) T increasing (decreasing) with system size 
at small (large) biases. 38 This renders this bias regime the most 
difficult for the determination of steady-state regimes. 

We emphasize the excellent agreement with the fRG 
results, 22 in particular, for bias voltages V > Tk- The differ- 
ences between tDMRG and fRG seen at larger biases for the 
U = t curve are due to the wide-band limit taken in fRG with 
a flat DOS (linear dispersion), compared to a semi-elliptical 
DOS in tDMRG calculations, as discussed above. We further 
compare to the 4th-order perturbation theory results (Ref. [13; 
dotted lines in Fig . O - While overall, the three methods agree 
for V > Tk and U/T = 2, the perturbation theory result lies 
slightly above both tDMRG and fRG for U/T > 4 and in the 
interval 0.2(7 < V < 0.8C7 (see Fig. |2i Note, though, that 
the differences between fRG and Ref. [lj are, at U/T = 4, 
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FIG. 4: (Color online) Mixed-valence regime: Current-voltage 
curves for U/T = 2 (U = t, symbols with dotted lines), U/T = 4 
[U = t/2, extrapolated data only, panels (b) and (d)], and several 
gate potentials: (a) V g = -17/4; (b) V g = 0; (c) V g = U/i; and (d) 
V g = U /2. Solid diamonds are extrapolations of finite-size data in 
l/N. Inset in (c): (n dot (r = 0)) vs. V s /U for U/T = 2 (circles) 
and 4 (squares). 



rather small and difficult to resolve numerically, but become 
more significant at U/T = 6 and 8. 

Therefore, we also consider the behavior for U/T = 8 and 
0.3 < V/U < 0.9, for which the accuracy of our data is 
reduced to S(J) T /V « 0.03Go, due to an additional oscilla- 
tory component in the transient behavior that becomes more 
prominent at large U /r4& Nevertheless, as Fig. |2jb) clearly 
shows, our tDMRG results are still in good agreement with 
fRG, while the corresponding 4th-order perturbation theory 49 
predicts too high values of (J) T /V. This constitutes an exam- 
ple where our numerical results serve to support the validity 
of one analytical approach, the fRG in this case, while unveil- 
ing that 4th-order perturbation theory^ 4 * seems to break down 
at U/T > 6 [see Fig.[2{b)]. Consequently, we may conclude 
that the double-minimum structure that Ref. [3 predicts to ex- 
ist in the differential conductance likely is an artifact of ap- 
proximations. This double minimum is precisely seen in the 
bias window in which the differences between our tDMRG 
and the 4th-order perturbation theory are the largest. 

As for the regime V < T K , note that at V/U = 0.1 and 
U/T = 2, the N — 96 result reaches about 99 % of the ex- 
pected perfect conductance. At larger U/T and V < 0.2U, 
the currents from < 128 are below the fRG result. Yet, for 
U/T = 4 and, e.g., V/U = 0.1 and 0.16, we have extrapo- 
lated the finite-size data for (J) T /V in l/N [see Fig. |2b)]. 
The agreement with fRG is quite reasonable. 



IV. RESULTS FOR THE MIXED- VALENCE REGIME 



We now turn our attention to the mixed-valence regime 
V g 7^ — U /2, concentrating on U/T = 2 and 4. Real-time cur- 
rents are shown in Fig.[3]for U/T = 2 (U = t) and V g = 0. 
In contrast to the V„ 



-U/2 case, significant oscillations 



with a period T oc \/V (Ref. |48I) in the current are observed. 
These are due to finite-size gaps and their amplitudes decay 
with 1/A. In principle, similar oscillations can be seen in the 
particle-hole symmetric case in the two local currents con- 
necting the dot to the left and right lead. As in computing 
(J) r , we average over these two local currents and as they 
oscillate with a phase shift of the oscillations cancel out 
in { J) T at V g = —U /2, but not away from that special point. 
For a detailed discussion of such oscillations in tDMRG cal- 
culations of real-time currents, see Ref. |4f| We exploit the 
existence of these oscillations to determine intervals to com- 
pute steady-state currents from, 4 2 averaging over a full period 
(see the arrows in Fig.[3]and the fit to the V = 1.2(7 data for 
an illustration). 

The average currents ( J) T are plotted vs. bias in Fig. [4] 
for U/T — 2 and 4. For U/T = 2, we include the results 
from system sizes of N = 32 and N = 64 sites to illustrate 
the emergent finite-size effects. To our advantage, we find 
that at large biases, even small chains of 32 sites yield well- 
converged results, as the example of V g = — U/4 and U/T = 
2 shows [see Fig. Ufa)]. Although we have not pursued this, it 
is reasonable to expect that at bias values of V/U > 2, exact 
diagonalization may produce useful results as well. 

For several values of the bias, we have extrapolated the av- 
erage currents in 1/A (shown as solid diamonds in Fig.|4]i. As 
compared to the V g = —U/2 case, finite-size effects are some- 
what more significant, as in our set-up, the average filling of 
the leads slightly exceeds hall-filling whenever V g ^ —U/2. 
The figure further includes the respective fRG results^ with a 
very good agreement with our extrapolated data. Slight devi- 
ations between fRG and tDMRG at intermediate bias voltages 
can be attributed to: (i) uncertainties in the extrapolation of 
the finite-size data and (ii) the curvature of the bands. Note 
that band effects already yield small differences in the linear 
conductance as obtained from tDMRG when compared to re- 
sults valid in the wide-band (linear dispersion) limit (see, e.g., 
Fig. 13 in Ref. ET 



V. SUMMARY 

In this work, we carried out large scale and accurate time- 
dependent simulations using the adaptive tDMRG method to 
study transport in the single-impurity Anderson model at large 
biases. By performing a finite-size analysis with chains of 
up to A = 128 sites, we were able to compute the current- 
voltage characteristics at the particle-hole symmetric point up 
to U/T = 8. Convergence with system size is very fast in 
the voltage regime Tk ^ V < 2t, while the resolution in 
the Kondo regime, i.e., V < Tk, is hampered by finite-size 
effects, due to the exponentially large Kondo correlations-^ 
(see Ref. [3^, though). Our results are in excellent agreement 
with those from fRGS for U/T < 8, lending strong sup- 
port to the validity of both approaches. Where converged, our 
data thus provide an unbiased benchmark from a quasi-exact 
method for analytical or approximate numerical approaches. 
Furthermore, we studied the mixed- valence regime by varying 
the gate potential at fixed U/T, presenting results for current- 
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voltage characteristics. 

Let us conclude by comparing our numerical approach, the 
tDMRG, to other numerical methods recently applied to non- 
equilibrium transport in the single-impurity model, namely, 
real-time QMC simulations, and the real-time NRG. A direct 
comparison of our results with the QMC data of Ref. Hg] (not 
shown here) reveals a reasonably good agreement at large bias 
voltages V > Tk- At small bias voltages, the QMC data are 
consistently below both tDMRG and fRG This can be traced 
back to the fact that, due to the sign problem, QMC simula- 
tions are limited to short time scales with a maximum time 
on the order of 3 /r^2 while in our case, we are able to reli- 
ably simulate times on the order of r < 6/r at large biases 
V ~ t and t < 12/r at small biases V ~ 0.2t (compare 
Fig. [TJ. We should emphasize, though, that the QMC simu- 
lations start from the dot being decoupled from the leads in 
the initial state, which is different from our initial condition 
and which may lead to a transient behavior different from our 
case. An advantage of the real-time QM C 29,30 is that it is de- 
signed for the thermodynamic limit, yet we have shown here 
that a finite-size analysis of tDMRG yields very accurate re- 
sults. Within both methods, the transient regime, i.e., the time 
window between the initial state and the steady-state regime 
can be studied4S 

Time-dependent NRG results were published for the dif- 
ferential conductance SG and large ratios of U/T (Ref. [l6h . 
One the one hand, one may generally expect time-dependent 
NRG to be well suited to capture the finite-bias behavior in the 
Kondo regime V < Tk, while on the other hand, the data pre- 
sented for SG at discrete bias voltages, seem to be somewhat 
noisy at large and intermediate biases (see Fig. 3c in Ref.[l6h. 
rendering a comparison with our data difficult. Real-time 



NRG calculates the steady-state density operator and then ex- 
tracts the current and the nonequilibrium spectral density from 
that information. The numerical challenges in obtaining accu- 
rate data for steady-state currents seem to be (a) keeping a 
sufficiently large number of states and (b) artefacts due to the 
logarithmic discretization. 16 Real-time NRG simulations can 
be performed at finite temperatures as wellr^ while methods 
for real-time DMRG simulations at finite temperatures are un- 
der developmen t 5 1 1 52 i 53 - 54 and may be available for the study 
of nonequilibrium currents in the near future. 

Our work establishes an important application of the adap- 
tive tDMRG to a topic of much interest to both experimen- 
talists and theoretical researchers. As the technique is highly 
flexible and can easily be adapted to impurity problems be- 
yond Eq. (fTJ such as interacting leads^ 2 - or complex geome- 
tries, we are convinced that this tool will become highly in- 
strumental in future studies of nonequilibrium transport in 
nanoscopic systems. 
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